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Abstract. The filamcntation instability of counterpropagating symmetric beams of 
electrons is examined with ID and 2D particle-in-cell (PIC) simulations, which are 
oriented orthogonally to the beam velocity vector. The beams are uniform, warm and 
their relative speed is mildly relativistic. The dynamics of the filaments is examined in 
2D and it is confirmed that their characteristic size increases linearly in time. Currents 
orthogonal to the beam velocity vector are driven through the magnetic and electric 
fields in the simulation plane. The fields are tied to the filament boundaries and 
the scale size of the flow-aligned and the perpendicular currents are thus equal. It 
is confirmed that the electrostatic and the magnetic forces arc equally important, 
when the filamcntation instability saturates in ID. Their balance is apparently the 
saturation mechanism of the filamcntation instability for our initial conditions. The 
electric force is relatively weaker but not negligible in the 2D simulation, where the 
electron temperature is set higher to reduce the computational cost. The magnetic 
pressure gradient is the principal source of the electrostatic field, when and after the 
instability saturates in the ID simulation and in the 2D simulation. 



PACS numbers: 52.40.Mj,62.35.Mw,52.65.Rr 
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1. Introduction 

The filamentation instability (FI) results in the growth and amplification of magnetic 
fields in astrophysical plasmas and in laser-generated plasmas. Magnetic fields grow, 
for example, due to the redistribution of currents in an initially current-free system of 
two counterstreaming electron beams [U [2] . This is the simplest plasma configuration 
that gives rise to the FI. The electronic FI grows faster than the competing two-stream 
and mixed mode instabilities, if the beams have a similar temperature and density and 
if their flow speed is at least mildly relativistic [HI HI [3]. A full classification of the 
competing modes in parameter space is given by Ref . j6j [7] . 

A wide range of previous numerical investigations of the FI exist, for example the 
pioneering PIC simulation studies in ID and in 2D systems [H [2]. Other important 
studies are the ID Vlasov [8] and the 2D PIC simulations [9] of counterstreaming 
electron beams, involving mobile ions. Relativistic beams and the impact of binary 
collisions on the FI have been investigated [10] . The impact of a guiding magnetic field 
on the FI has been studied [HI [121 03] as well as the combination of filamentation and 
Weibel instabilities [14]. Equilibrium conditions of the beam-plasma system have been 
addressed pj)]. The statistical distribution of filament sizes has been the focus of ID 
and 2D PIC simulation studies [16-19], while the astrophysical relevance of the FI driven 
by leptonic beams has been assessed in the Refs. [3, 20-26]. 

Here we revise the special case of the FI, which is driven by counterstreaming 
electron beams that are equally dense and moderately warm. This FI is electromagnetic 
during its linear growth phase [H [TJ. We model this FI with particle-in-cell (PIC) 
simulations. The growth and the saturation of the FI can be described as follows. 
Initially, the currents of the counterpropagating electron beams cancel each other. The 
noise inherent to PIC simulations results in fluctuating magnetic fields. A magnetic field 
fluctuation gives rise to a small separation of the electrons of both beams and, thus, to 
a net current. This net current amplifies, in turn, the magnetic field. The latter then 
grows exponentially until the FI saturates [1]. 

Simulations have evidenced the nonlinear growth of electrostatic fields by the FI 
and their importance has been pointed out [8j [10] . The comparison of the electric and 
magnetic field profiles [131 EH] for the case of symmetric beams, when no electrostatic 
field can grow in the linear phase of the FI, suggested, that the source mechanism is 
the magnetic pressure gradient force (MPGF). It induces through the acceleration of 
electrons a current, which drives the electrostatic fields. The fields oscillate around 
an equilibrium amplitude set by the MPGF [27] after the FI has saturated, unless 
positrons are present [26]. This current and the electrostatic fields are not affected by 
the introduction of a spatially uniform flow-aligned magnetic field due to its vanishing 
contribution to the MPGF [13]. A flow-aligned uniform magnetic field apparently only 
reduces the linear growth rate of the FI or it suppresses it [TTl [T3j . 

This paper is structured as follows. Section 2 outlines the simulation parameters. 
Section 3 revises a study [27] of the special case of symmetric cool electron beams. 
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The link between the saturated electrostatic field and the MPGF is illustrated and it 
is shown that the FI saturates due to an approximate cancelling of the magnetic and 
electrostatic forces. It is shown quantitatively, rather than qualitatively p3], for the 
first time that the saturated electric field in the simulation plane equals that expected 
from the MPGF in a 2D simulation. Section 5 is the summary. 



2. The PIC simulation method and the initial conditions 

The PIC method can, in principle, model all plasma processes in a collisionless plasma. 
It approximates the plasma phase space distribution by an ensemble of computational 
particles (CPs), each of which has a position x cp and velocity v cp . Their charge to mass 
ratio equals here the —e/m e of an electron. The ions form in the simulations discussed 
here an immobile charge background, which cancels the net electron charge. 

A PIC code solves the discretized Maxwell equations for the fields and the Lorentz 
equation for each CP and it interpolates the quantities defined on the grid to the 
positions of the particles and vice versa [28]. Our numerical scheme is outlined in 
Ref. [29]. The physical quantities are normalized as follows. The plasma frequency 

1 /2 

Up = (2n e e 2 /m e eo) is obtained from the summed density of the two electron beams 
we model. Each beam has initially the spatially uniform density n e . The skin depth 
A e = c/u p . The quantities in physical units denoted by the subscript p are obtained from 
the normalized ones by substituting E p = uj p cm e E / e, B p = u p m e B/e, J p = 2en e cJ, 
p p = 2en e p, x p = \ e x, t p = t/u p , v p = vc and p p = m e cp. We also normalize fl = uj /uj p 
and k = k p c/uj p , where uj and k p are given in physical units. The equations are 

VxE = -—, VxB = — + J,V-B = 0,V-E = p, (1) 



dp C p 



q cp (E[x cp \ + v cp x B[x cp }) . (2) 



dt 

Initially E = and B = 0. The electron beams have each a speed modulus vj, = 0.3c 
and they move in opposite directions along z. The thermal speed v e = (kBT e /m e ) ' 5 
is set to Vb = 9v e (18t> e ) in the 2D (ID) simulation. The electron temperature in 
the 2D simulation is higher to reduce the computational cost by the increased Debye 
length, which determines both the grid cell size and the maximum time step that is 
possible. We resolve the x direction (ID) or the xy plane (2D), by which we isolate the 
FI with its k _L z. The 2D simulation uses 1500 x 1500 grid cells to resolve the domain 
L x x L y = 90A e x 90A e . Each electron beam is represented by 144 CPs per cell. The ID 
simulation resolves Li = 0.89A e by 500 grid cells and each electron beam by 1.21 x 10 5 
CPs per cell like in Ref. [27]. The boundary conditions are periodic. 



3. Simulation results 



The FI results in the separation of the currents C z (x,y) along z of both beams. 
We examine it and the current component in the (x,y)-plane, which we denote with 
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(a) Position X, Current C (b) Position X, Current C 



Figure 1. (Colour online) The current in a section of the 2D box at t = 33, normalized 
to the initial current of one beam: (a) displays the flow-aligned current C z and (b) the 
modulus of the current C xy = C x +iC y in the simulation plane. No obvious connection 
exists between the beam-aligned current and the perpendicular current. 

C xy (x, y) = C x (x, y) + iC y (x, y). Both currents at the time t = 33 in the 2D simulation 
are displayed in Fig. (TJ when the FI has just saturated. Both beams and their 
contributions to C z (x, y) have been separated into domains. The peak modulus is about 
twice the mean current of one beam. The structures in \C xy \, which have a significant 
strength, show in some cases a correlation with those in C z . The evolution of C z (x, y) is 
animated in time in the movie 1 and that of \C xy (x, y)\ in the movie 2. Movie 1 shows 
the initial growth of stationary filaments. These start to merge and deform during the 
nonlinear phase of the FI. The filament dynamics slows down in time, as the filament 
size increases relative to the boundary speed. Movie 2 demonstrates, how structures in 
\C xy \ come and go. They are damped and must thus be driven by C z . 

The C z (x,y,t) and the C xy (x,y,t) are Fourier transformed to C z (k x , k y ,t) and 
C xy {k x , k y ,t) . The P z [k x ,k y ,t) = \C z {k x , k y , t)\ and P xy {k x ,k y ,t) = \C xy {k x , k y ,t)\ 
are calculated. We transform (k x ,k y ) — > {k cos a, k sin a) and integrate the power over 
a in the k-plane to give the P z (k, t) and P xy (k, t) shown in Fig. [21 The P z (k, t) is stronger 
than P xy (k,t), as expected. The k-interval, in which the power peaks in both current 
components, shifts in time like k oc t -1 . The characteristic filament size oc k~ x thus 
increases linearly with the time. It is evident that the peak power of both spectra is 
located in the same k-interval at any fixed time, after the FI has saturated. The power 
P X y(k,t) maintains its slope at high k, while that of P z (k,t) is broadening in time. The 
shift of the boundary of P z (k, t) at low k to larger values of k at late times might be a 
finite box effect, as we notice the discreteness of k. 

Figure [3] displays B x , B y and the normalized magnetic pressure Pb = + By)/2 
at t = 235. The magnetic field vanishes within the filaments and it is strong at their 
boundaries, as it is demonstrated by Pb(x, y). The magnetic reconnection points, e.g. 
at x = 12 and y = 20, and the locations (x = 20, y = 20), where reconnection is about to 
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Figure 2. (Colour online) The spatial power spectrum of the currents in the 2D box, 
integrated over the azimuth in (k x , k y ) as a function of time: (a) corresponds to the 
P z (k, t) and (b) to P xy (k, t). Both power spectra are normalized to the maximum in (a) 
and the colour scale is 10-logarithmic. The same two curves k cx t -1 are ovcrplottcd. 




X, Timet = 235 (c) X, Time t = 235 



Figure 3. (Colour online) The magnetic fields at t = 235 in the 2D simulation: (a) 
displays the B x in a subsection of the simulation box and (b) the B y in the same 
interval. The magnetic pressure Pb = (B x + By)/ 2 in the full box is shown in (c). 

take place, demonstrate the merging of filaments. The Pb(x, y) reveals strong gradients, 
which should result in a significant MPGF. A ID PIC simulation provides more insight 
into the relevance of the MPGF. This can be exemplified with the help of the fluid 
momentum equation [30], which we will reduce for this purpose to one dimension x. 
Each species with the index s is described by a fluid with the density n s , the mean 
speed v s and the pressure tensor P s and the equation of motion in SI units is 

q n 

d t {n s v s ) + V • (n s v s v s ) = V • P s + — (E + v s x B) . (3) 

m e m s 

Ampere's law is used to rewrite J s x B with J s = q s n s v s . We obtain the equation for 
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the electron species s — 1 and s = 2, which we model. 

d t (n s v s ) + V • (n s v s v s ) = -— • P s - ^£ - J-^ - — B x «9 t £ (4) 

The contribution by the gradient of the magnetic stress tensor [1q 1 (BB) is omitted. 
Its contribution is zero if we go to a ID geometry, because the gradients along y and z 
vanish and because B x (x,t) = 0. The thermal pressure tensor is diagonal for our initial 
conditions. The FI results in the ID simulation in the initial growth of B y and E z . The 
electrostatic E x grows nonlinearly. All other field components remain at noise levels. 
The component of Eq. H] along x simplifies to 

a i \ . j / 2 \ o n sk B T s en s B v d x B y e 

d t {n s v StX ) + d x {n s v sx ) = -d x E x ■ B y d t E z \b) 

m e m e m e p m e 

The electric and magnetic fields are computed self-consistently by the PIC simulation 
using the total charge and current density, which are obtained as p = p\ + p 2 + p% 
and J = J\ + J2 from the electron species 1 and 2 and from the background ions. The 
uniform and constant charge density p\ of the immobile ions cancels initially the electron 
charge and the ions do not provide a current. The d x (n s Vg X ) and the thermal pressure 
gradient can be neglected during the linear growth phase, when only B y and E z grow. 
We retain d t {n s v s ^ x ) oc B y d x B y . The displacement current can probably be neglected, 
because of B y /E z « 100 [27| and because E z grows smoothly and aperiodically. A 
current grows along x in response to the MPGF, which drives an electrostatic field 
through d t E x = —J x /e or dfE x oc —B y d x B y , because V x B = in a ID simulation. 

The right hand side of Eq. [5]is zero if en s E x = — B y d x B y / p for the species s, if the 
thermal pressure gradient and the term with the displacement current are negligible. 
This may provide a condition for the saturation of the FI. The aperiodic growth of E x 
may be prescribed by dfE x oc —B y d x B y and the spatial amplitude of E x should in this 
case be proportional to that of the MPGF. An exact cancelling of the term oc E x and 
the MPGF in Eq. [5]is not possible, if n s varies as a function of x. The equations of both 
electron fluids are summed up to get a condition (n\ + n 2 )E s = —2B y d x B y / p for the 
saturation electric field E$ and we find that 711(2;) +77,2(2;) ~ 2n e for this particular case 
study [27]. The normalization (section 2) to E$ = —2B y d x B y facilitates the comparison 
of this feasible saturation condition with what we observe in the simulation. The factor 
2 arises from the normalization to 2n e in section 2. We define Eb = Eg/2. 

Figure H] compares the fields computed by the ID PIC simulation. Only one wave 
period of B y (x,t) is resolved by the box length L\. The B y saturates at t ~ 56 and it 
remains approximately stationary thereafter. The saturation of B y is accompanied by 
the growth of E x . The E x has twice the wavenumber of B y and it oscillates around 
a background electric field, which is stationary in space. The latter has the same 
wavenumber and amplitude as Eb- The Eb and E x show correlations, e.g. at t — 63. We 
demonstrate the quantitative match of 2E B (x, t = 56) and E x (x, t = 56), when the latter 
reaches its peak amplitude. This accurate match confirms the dominance of the MPGF 
over the displacement current term at this time, when the second and third term on the 
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(a) Position X ( b ) Position X ;c ) Position X £e) Potion x 



Figure 4. (Colour online) The field data from the ID simulation: (a) and (b) show 
B y (x,t) and E x (x,t), respectively. The is displayed in (c). The E x (x,t = 56) 

(solid) and 2EB(x,t = 56) (dashed) are compared in (d), while (e) compares the 
time-averaged Exa(x) (solid)and Eba(x) (dashed). 

right hand side of Eq. [5] practically cancel each other. A good agreement is also achieved 
between Eba{ x ) = {h — ^i)" 1 / Eb(x, t) dt and Exa{%) = (£2 ~ ^i)" 1 / E x (x, t) dt, which 
have both been integrated from t\ = 60 to £2 = 130. The E x (x,t) oscillates with 
the spatial amplitude Eba{%) around the stationary background field Eba{x). They 
add up to E x (x,t = 56) = 2Eba(x) and they cancel each other at the times to, when 
E x (x, to) = in Fig. H](b). Ions will only weakly react to the high-frequency oscillation of 
E x , but the stationary Eba will accelerate them [8]. The ion current and its nonuniform 
charge will eventually modify the fields and the force balance. 

Figure [5] compares the impact of E x and of B y on the saturation of the FI in the 
ID simulation and on the phase space distribution of the beam moving with v = +u& for 
the two times t = 45, 56. The phase space distributions demonstrate that all electrons 
have a v z ~ v b and a \v x \ Vj,. Movie 3 demonstrates that this is true throughout the 
simulation by animating the phase space densities log 10 f(x, v z , t) and log 10 f(x, v x , t). 
This weak heating may be the reason, why the thermal pressure gradient has no obvious 
influence on the saturation of the FI. Movie 3 reveals the rearrangement of the electrons 
of this beam into a filament and the vortex formation in f(x,v x ). 

The normalized electric force on an electron is — E x . The normalized magnetic 
deflection force working on an electron of the beam with Vb > is ~ —VbB y . The 
magnetic force exceeds by far the electric one at t = 45. Both forces are comparable at 
t = 56, when the FI saturates, which serves as a further illustration for the saturation 
condition obtained from Eq. [5j The magnetic deflection force is responsible for the 
trapping of the filament electrons. The E x is repulsive at the centre of the filament 
and counteracts the magnetic trapping. It is attractive at longer distances, facilitating 
the filament overlap [27] . The E x limits the peak density of the filament and, thus, the 
current it can carry [26] . The combined action of both forces is to confine the electrons 
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Figure 5. (Colour online) Projections f(x, v z ) and f(x, v x ) of the electron phase space 
distributions and plots of the fields at t = 45 (left column) and t = 56 (right column). 
The colour scale of f{x,v z ) in (a,d) and f(x,v x ) (b,e) is linear and in units of CPs. 
The curves VbB y (solid) and E x (dashed) are plotted in (c,c). 




(a) X, Time t = 235 (b) X, Time t = 235 

Figure 6. (Colour online) The magnitude of Eb is shown for a subsection of the box 
at t = 235 in (a). The electric field modulus \E X +iE y \ is displayed in (b) for the same 
time and box interval. Both distributions have been computed by the 2D simulation. 



into a filament in f(x,v z ) and a vortex in f(x,v x ). The periodic oscillations of these 
electrons and their J x in the potential (Movie 3) give rise to the oscillatory E x . 

Figure [6] compares the magnitude of Eb = — (d x B 2 ,d y B 2 ) with the modulus of 
the complex electric field E x + %E y in the 2D simulation at t — 235. We select a late 
time, because then the filament dynamics is not so fast (movie 1) compared to the 
oscillation frequency of the electric field (Fig. @|. The boundaries are also quasi-planar, 
by which they become locally one-dimensional. This should reduce the importance of 
the magnetic tension force relative to the MPGF. The \Eb\ and \E X + iE y \ are spatially 
correlated, they have the same magnitude and they reveal a split in the centre of the 
band. This split occurs in the ID simulation, when d x B y = 0. It is thus likely, that 
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the electric field in the simulation plane is driven practically exclusively by the MPGF. 
Both fields are out of phase at some places, e.g. at (x,y) = (10, 10), where \E B \ does 
not have a maximum. This may be explained with the interplay of J with E, causing 
oscillations of the electric field. The peak modulus of \Eb\ ~ 0.04, which is about half 
the peak value of the magnetic deflection forces f&l-Rrl or Vb\B y \ in Fig. [3j 

4. Discussion 

We have investigated the FI with ID and 2D PIC simulations, where the beam velocity 
vector is orthogonal to the simulation direction or plane. The uniform, cool and 
symmetric electron beams have a mildly relativistic relative speed. The 2D simulation 
has provided insight into the interplay of the filaments, e.g. their merging, during the 
nonlinear phase of the FI. The filaments formed by the beam-aligned current couple to 
the perpendicular currents through the electromagnetic fields. The fields are confined to 
the filament boundaries and the characteristic size of the structures in the flow-aligned 
and the perpendicular current are thus equal. This characteristic size increases linearly 
with the time. The power spectrum of the damped perpendicular currents showed for all 
times after the FI has saturated a constant slope at large wavenumbers, which follows 
approximately a power law, albeit over a narrow range of wavenumbers [18] . 

The ID PIC simulation using the same setup as in Ref. [27] revealed a feasible 
source mechanism for the electrostatic fields, which develop during the nonlinear 
phase of the filamentation instability. The MPGF accelerates the electrons and a 
current develops, which couples through Amperes law into the electrostatic field. This 
mechanism is suppressed by selecting beams of equally dense electrons and positrons, 
because here the currents of both species cancel |26j . The growth of the electrostatic 
fields is not affected by a spatially uniform and flow-aligned magnetic field, because its 
MPGF contribution vanishes [13]. It can reduce the growth rate of the FI though. 

The current and the electric field oscillate in the ID simulation after the FI has 
saturated. The spatial profile of the time-averaged electrostatic field amplitude is that 
expected from the MPGF. The electric field amplitude oscillates between the initial 
value E x (x, t = 0) = and twice its mean amplitude. The nonlinear terms in a ID fluid 
equation due to the electrostatic field and the MPGF cancel approximately, when the 
electrostatic field peaks. This may suggest that, for the special case of symmetric and 
nonrelativistic electron beams considered here, the FI saturates due to the balancing 
of the electrostatic and magnetic forces. We demonstrated for the first time that the 
MPGF is also responsible for the electric field in the simulation plane of the 2D PIC 
simulation. The electric field amplitude in the 2D simulation is weaker than that in the 
ID simulation, but it is not negligible. This might be a consequence of the increased 
electron temperature that reduced the computational cost of the 2D simulation. The 
magnetic energy density we obtain is about 5% (10%) of the total energy in the ID (2D) 
simulation. The magnetic field amplitude for a plasma frequency « 10 kHz, which is 
representative for the solar wind, would be about 20 nT. 
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Our findings will be relevant for (almost) symmetric electron beams. Selecting 
asymmetric beams facilitates plasma equilibria [15], rather than a time-dependent 
evolution. The FI will then also result in the generation of electrostatic fields during 
its linear growth phase [I] and the mixed mode instability will outrun the FI. Highly 
relativistic beams will probably break the simple relation between the MPGF, the 
current and the electrostatic field. The thermal pressure gradient and the magnetic 
tension may not be negligible for other initial conditions and during the initial nonlinear 
stage, when the filament boundary curvature is higher. The competition of the FI with 
the mixed mode instability and the electrostatic two-stream instability, which saturates 
by forming phase space holes, as well as the flux tube bending, will introduce nonlinear 
effects that can only be addressed with large-scale 3D PIC and Vlasov simulations. 
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